#include <iostream>
#include<sys/time.h>
#include<omp.h>
#include<arm_neon.h>
#define N 500
#define C 1
#define P 7
using namespace std;
float A[N][N];
void display() {
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            cout<<A[i][j]<<" ";
        }
        cout<<endl;
    }
}//check outcome
void init() {
    //cout<<"Initializing ...\n";
    for(int i = 0; i < N; i++) {
        A[i][i] = 1.0;
    }
    for(int r = 0; r < 5 * N; r++) {
        int i1 = rand()%N;
        int i2 = rand()%N;
        float rate = rand()%10/10.0;;
        if(i1 != i2) {
            for(int j = 0; j < N; j++) {
                A[i1][j] += rate*A[i2][j];
            }
        }
    }
    //cout<<"Finished. \n";
}
//平凡算法
void normal() {
    struct timeval h, t;
    float time = 0.0;
    //Normal
    cout<<"Normal:(N = "<<N<<")\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            for(int i = k + 1; i < N; i++) {
                for(int j = k + 1; j < N; j++) {
                    A[i][j] = A[i][j] - A[i][k]*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
//openMP
void om() {
    struct timeval h, t;
    float time = 0.0;
    //Normal
    cout<<"ROW:(N = "<<N<<", P = 3)\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        int i,j,k;
        float tmp;
        #pragma omp parallel num_threads(3),private(i, j, k, tmp)
        for(k = 0; k < N; k++) {
            #pragma omp master
            {
                tmp = A[k][k];
                for(j = k + 1; j < N; j++) {
                    A[k][j] = A[k][j]/tmp;
                }
                A[k][k] = 1.0;
            }
            #pragma omp for schedule(guided, 5)
            for(i = k + 1; i < N; i++) {
                tmp = A[i][k];
                for(j = k + 1; j < N; j++) {
                    A[i][j] = A[i][j] - tmp*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
void om1() {
    struct timeval h, t;
    float time = 0.0;
    //Normal
    cout<<"ROW:(N = "<<N<<", P = 7)\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        int i,j,k;
        float tmp;
        #pragma omp parallel num_threads(7),private(i, j, k, tmp)
        for(k = 0; k < N; k++) {
            #pragma omp master
            {
                tmp = A[k][k];
                for(j = k + 1; j < N; j++) {
                    A[k][j] = A[k][j]/tmp;
                }
                A[k][k] = 1.0;
            }
            #pragma omp for schedule(guided, 5)
            for(i = k + 1; i < N; i++) {
                float tmp = A[i][k];
                for(j = k + 1; j < N; j++) {
                    A[i][j] = A[i][j] - tmp*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
void om2() {
    struct timeval h, t;
    float time = 0.0;
    //Normal
    cout<<"COL:(N = "<<N<<", P = 7)\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        int i,j,k;
        float tmp;
        #pragma omp parallel num_threads(7),private(i, j, k, tmp)
        for(k = 0; k < N; k++) {
            #pragma omp master
            {
                tmp = A[k][k];
                for(j = k + 1; j < N; j++) {
                    A[k][j] = A[k][j]/tmp;
                }
                A[k][k] = 1.0;
            }
            #pragma omp for schedule(guided)
            for(j = k + 1; j < N; j++) {
                float tmp = A[k][j];
                for(int i = k + 1; i < N; i++) {
                    A[i][j] = A[i][j] - tmp*A[i][k];
                }
            }
            for(i = k + 1; i < N; i++) {
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}
void om3() {
    struct timeval h, t;
    float time = 0.0;
    //Normal
    cout<<"ROW+NEON:(N = "<<N<<", P = 7)\n";
    for(int r = 0; r < C; r++) {
        init();
        gettimeofday(&h, NULL);
        int i,j,k;
        float tmp;
        #pragma omp parallel num_threads(7), private(i, j, k, tmp)
        for(k = 0; k < N; k++) {
            #pragma omp master
            {
                tmp = A[k][k];
                float* akk = new float[4];
                *akk = *(akk+1) = *(akk+2) = *(akk+3) = tmp;
                float32x4_t vt = vld1q_f32(akk);//将四个单精度浮点数从内存加载到向量寄存器
                delete[] akk;
                for(j = k + 1; j + 3 < N; j += 4) {
                    float32x4_t va = vld1q_f32(&A[k][j]);
                    va = vdivq_f32(va, vt);//A[k][j] = A[k][j]/A[k][k];
                    vst1q_f32(&A[k][j], va);//储存到内存
                }
                for(j = N - N%4; j < N; j++) {
                    A[k][j] = A[k][j]/tmp;
                }//结尾几个元素串行计算
                A[k][k] = 1.0;
            }
            #pragma omp for schedule(guided, 5)
            for(i = k + 1; i < N; i++) {
                tmp = A[i][k];
                float* aik = new float[4];
                *aik = *(aik+1) = *(aik+2) = *(aik+3) = tmp;
                float32x4_t vaik = vld1q_f32(aik);//将四个单精度浮点数从内存加载到向量寄存器
                delete[] aik;
                for(j = k + 1; j + 3 < N; j += 4) {
                    float32x4_t vaij = vld1q_f32(&A[i][j]);
                    float32x4_t vakj = vld1q_f32(&A[k][j]);
                    float32x4_t vx = vmulq_f32(vaik, vakj);
                    vaij = vsubq_f32(vaij, vx);//A[i][j] = A[i][j] - A[i][k]*A[k][j];
                    vst1q_f32(&A[i][j], vaij);//存储到内存
                }
                for(j = N - N%4; j < N; j++) {
                    A[i][j] = A[i][j] - tmp*A[k][j];
                }//结尾几个元素串行计算
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME: "<<time/C<<" ms"<<endl;
}

int main()
{
    

        normal();

        om();

        om1();

        om2();

        om3();

}